#ifndef __QGCM_PROFILE_MAPPERS__
#define __QGCM_PROFILE_MAPPERS__

#include "forward_classes.h"
#include "datatypes.h"
#include "readers.h"
#include "static.h"

/* -------------------------------------------------------------------------------- 


 */

class custom_func_container {
	model_class & mo;
	ParserX * cval; // universal parser rule 

	const string OPname;  

	const int rank_fields; // number of fields
	const int rank_kspace; // number of k space variables
	const int rank_amodes; // number of angular modes

	vector< pair<string,Value*> > fields;
	vector< pair<string,Value*> > kvariables;
	vector< pair<string,Value*> > angmodevariables;

	public:
	custom_func_container(const string &, const int, const int, const int);
	~custom_func_container();

	const IValue & eval();
	// include all overloaded pairs

	bool isdefined() const;

	void assignkspacenames( const vector<string> & );
	void assign_angspace_names( const vector<string> & );
	void assign_field_names( const vector<string> & );

	void assignkspace(REAL);
	void assignkspace(REAL, REAL);
	void assignkspace(REAL, REAL, REAL);

	void assignnkspace(INDEX, INDEX);

	void assignindices(const rank3 &, const rank3 &);
	void assignindices(const rank3 &);
	void assignindices(const rank2 &);

	void assignexpr(const token_list & expr);

	const string & name() const;

	int getnum_ang_modes() const;
	int getnum_kspace_variables() const;
	int getnum_fields() const;

	private:
	custom_func_container();
	custom_func_container(const custom_func_container &);
	custom_func_container & operator=(const custom_func_container &);
};
// --------------------------------------------------------------------------------
class get_custom_func_container {
	custom_func_container * wfunc;
	static custom_func_container asym_kpq_cp;
	static custom_func_container sym_kpq_cp;
	static custom_func_container kpq_cp;
	static custom_func_container R_cp;
	static custom_func_container S_cp;
	static custom_func_container customAhat;
	static custom_func_container customQ;
	static custom_func_container customP;
	public:
	get_custom_func_container(int); // get by index
	get_custom_func_container(const string &); // get by name
	custom_func_container & operator*(); 
	bool isvalid() const;
	bool operator!() const;
};
// --------------------------------------------------------------------------------
// custom container global functions
void custom_symcomputeKPQintegrand(REAL & rvalue, const REAL & k, const REAL & p, const REAL & q);
qgcomplex custom_asymcomputeKPQintegrand(const REAL & k, const REAL & p, const REAL & q, const INDEX & np, const INDEX & mk);
REAL Rcustom(const model_class & mo, REAL k2, REAL p2, REAL q2, size_t i1, size_t i2);
REAL Scustom(const model_class & mo, REAL k2, REAL p2, REAL q2, size_t i1, size_t i2);
REAL customQ(const model_class &, const REAL & k, const rank2 & field_indices);
REAL customAhatfunc(const model_class & mo, const REAL & k, const REAL & p, const REAL & q, const rank3 & ipos);
bool customP(const rank3 & ipos);

/* -------------------------------------------------------------------------------- 
   Below we define the mirror_array class that satisfies the property
   data[-i] = data[i]
   data[0 ] = data[ 0];
   Nothing special. It is a wrapper for the std::vector< .. > class
   and essentially redefines the operator[] to achieve this goal.
 */
// --------------------------------------------------------------- //
class mirror_array {
	vector<INDEX> * data;
	public:
	mirror_array(const vector<INDEX> & array);
	mirror_array();
	~mirror_array();
	mirror_array & operator=(const mirror_array &);
	mirror_array & operator=(const vector<int> &);
	size_t size() const;
	INDEX operator[](INDEX index) const;
	const INDEX & back() const;
};

/* ------------------------------------------------------------------------- */
// Function below generates the M nonlinear operator matrix. This is not
// used in the model anywhere, but is useful for diagnosing the accuracy
// of Q and Qinv
/* ------------------------------------------------------------------------- */


// ......................................................................... //
class computestaticpart {
	const model_class * mo;

	vector<rank3> Afieldindexmap;
	static_container * Rstaticweights;
	static_container * Sstaticweights;
	map_container<REAL, 6> * Ahatmap;
	map_container<REAL, 3> * Qmap;
	map_container<REAL, 3> * QImap;


	void ( * symkpqintegrand ) ( REAL &, const REAL &, const REAL &, const REAL & );
	qgcomplex ( * asymkpqintegrand ) ( const REAL &, const REAL &, const REAL &, const INDEX &, const INDEX & );

	REAL ( * func_compute_Q ) ( const model_class & mo, const REAL &, const rank2 & );
	REAL ( * func_compute_A ) ( const model_class & mo, const REAL &, const REAL &, const REAL &, const rank3 & );  
	//custom approach to computingA
	bool ( * func_compute_P  ) ( const rank3 & );

	const REAL R_coeff;
	const REAL S_coeff;
	const vector<int>  nresols;
	const vector<REAL> kdomain;
	const vector<REAL> dkdomain;
	const vector<REAL> kdomain_centered;
	const vector<REAL> kspacing;
	const bool ignoreslot;
	const int NUMPROCESSORS;
	const int numfields;
	const bool ignorebasetriangle;
	const REAL kbasefactor; 
	const REAL THRESHOLD;
	const REAL dk_width_adjuster;
	const mirror_array nkdomain;
	public: 
	computestaticpart(const model_class & model,
			static_container * & Rs,
			static_container * & Ss,
			map_container<REAL, 6> * & Amap,
			map_container<REAL, 3> * & Qmap,
			map_container<REAL, 3> * & QImap,
			void ( * kpqintegrand ) ( REAL &, const REAL &, const REAL &, const REAL & ),
			complex <REAL> ( * asymkpqintegrand )( const REAL & , const REAL &, const REAL &, const INDEX &, const INDEX &),
			bool ( * fcompute_P ) ( const rank3 & ),
			REAL ( * fcompute_Q ) ( const model_class & mo, const REAL &, const rank2 & ),
			REAL ( * fcompute_A ) ( const model_class & mo, const REAL &, const REAL &, const REAL &, const rank3 &),
			const REAL R_coeff,
			const REAL S_coeff
			);
	private:
	computestaticpart(); // turn off default constructor 

	REAL computeA( const INDEX &, const INDEX &, const INDEX &, const rank3 & ) const;

	// now for important functions
	void initialize();  // init all member data
	void mapstaticpart(); // actually computation of the static part

	void mapAhat(REAL * & functions);
	bool mapQ();

	void mapsympart(const size_t, REAL *, rank3 *) const; // wrapper func for fine_mesh_integrate that avoids double work 
	bool preassignsymvalue(REAL &, const REAL *, const rank3 *, const rank3 &, const size_t &) const;
	int findsymvalue( const rank3 *, const rank3 &, const size_t &) const;
	static rank3 sortrank3(const rank3 &);

	REAL fine_mesh_integrate(const INDEX ki, const INDEX pi, const INDEX qi, const float *) const;
	inline bool isinsideslot(const REAL, const REAL, const REAL) const;
	inline bool isinsideregularizedslot(const REAL, const REAL, const REAL) const;
#ifdef DEBUG
	static void initialize_symmetric_data_array(const size_t, REAL *, rank3 *);
	//#define isinsideslot(k,p,q) (( q < k && p <= k + q && p >= k - q ) || ( q >= k && p <= k + q && p >= -k + q))
#endif
};


/* ------------------------------------------------------------------------- */
#endif

